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Diffusion of transcription factors can drastically enhance the noise in gene expression 
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We study by simulation the effect of tlie diffusive motion of repressor molecules on the noise in 
mRNA and protein levels in the case of a repressed gene. We find that spatial fluctuations due to 
diffusion can drastically enhance the noise in gene expression. For a fixed repressor strength, the 
noise due to diffusion can be minimized by increasing the number of repressors or by decreasing the 
rate of the open complex formation. We also show that the effect of spatial fluctuations can be well 
described by a two-step kinetic scheme, where formation of an encounter complex by diffusion and 
the subsequent association reaction are treated separately. Our results also emphasize that power 
spectra are a highly useful tool for studying the propagation of noise through the different stages 
of gene expression. 
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I. INTRODUCTION 



Cells process information from the outside and regu- 
late their internal state by means of proteins and DNA 
that chemically and physically interact with one another. 
These biochemical networks are often highly stochastic, 
because in living cells the reactants often occur in small 
numbers. This is particularly important in gene expres- 
sion m, 1^ yf , where transcription factors are frequently 
present in copy numbers as low as tens of molecules per 
cell. While it is generally believed that biochemical noise 
can be detrimental to cell function |j], it is increasingly 
becoming recognized that noise can also be beneficial to 
the organism 5] . Understanding noise in gene expression 
is thus important for understanding cell function, and 
this observation has recently stimulated much theoreti- 
cal and experimental work in this direction J4| Q . How- 
ever, the theoretical analyses usually employ the zero- 
dimensional chemical master equation 0, Q]. This ap- 
proach takes into account the discrete character of the 
reactants and the probabilistic nature of chemical reac- 
tions. It does assume, however, that the cell is a 'well- 
stirred' reactor, in which the particles are uniformly dis- 
tributed in space at all times; the reaction rates only de- 
pend upon the global concentrations of the reactants and 
not upon the spatial positions of the reactant molecules. 
Yet, in order to react, reactants first have to move to- 
wards one another. They do so by diffusion, or, in the 
case of eukaryotes, by a combination of diffusion and ac- 
tive transport. Both processes are stochastic in nature 
and this could contribute to the noise in the network. 
Here, we study by computer simulation the expression of 
a single gene that is under the control of a repressor R in 
a spatially-resolved model. We find that at low repressor 
concentration, i. e. [R] < 50nM, the noise in gene expres- 
sion is dominated by the noise arising from the diffusive 



motion of the repressor molecules. Our results thus show 
that spatial fluctuations of the reactants can be an im- 
portant source of noise in biochemical networks. 

Our analysis reveals that in gene expression significant 
fluctuations occur on both small and large length and 
time scales. As expected from earlier work |3, 0, U^, 
the fluctuations on long time scales are predominantly 
due to protein degradation; we assume that proteins are 
degraded by dilution, which means that the relaxation 
rate of this process is on the order of an hour. Our re- 
sults, however, also elucidate an important process on 
much smaller length and time scales. It is associated 
with the competition between repressor and RNA poly- 
merase (RNAP) for binding to the promoter. When a 
repressor molecule dissociates from the DNA, it can re- 
bind very rapidly, i.e. on a time scale of microseconds, or 
less. This rebinding time is so short that when a repres- 
sor molecule has just dissociated, the probability that a 
RNAP will bind before the repressor molecule rebinds, 
is very small. As a result, a repressor molecule will on 
average rebind many times, before it eventually diffuses 
away from the promoter and a RNAP molecule, or an- 
other repressor molecule, can bind to the promoter. This 
process of rapid rebindings decreases the effective dissoci- 
ation rate, and this increases the noise in gene expression. 

Clearly, fluctuations in gene expression span orders of 
magnitude in length and time scales. This means that 
the simulation technique should be sufficiently detailed 
to resolve the events at small length and time scales, 
yet also efficient enough to access the long length and 
time scales. Recently, several simulation techniques have 
been developed for the stochastic modeling of reaction- 
diffusion systems [lll[l2- These techniques, however, do 
not satisfy both criteria: they either describe the system 
in a coarse-grained way, i. e. on the level of local concen- 
trations rather than single particles ^^ [l3 , or are too 
slow to accurately model the dynamics on the long time 



scales [13] ■ Our simulations have been made possible via 
the use of our recently developed Green's Function Re- 
action Dynamics (GFRD) algorithm HI Ell • GFRD is 
an event driven algorithm that uses Green's functions to 
combine in one step the propagation of the particles in 
space with the reactions between them. The event-driven 
nature of the algorithm makes it particularly useful for 
problems, such as gene expression, in which the events 
are distributed over a wide range of length and time 
scales: the algorithm takes small steps when the reac- 
tants are close to each other - such as when a repressor 
molecule has just dissociated from the DNA - while it 
takes large jumps in time and space when the molecules 
are far apart from each other - like when the repressor 
molecule has eventually diffused away from the promoter. 
The event-driven nature of GFRD makes it orders of 
magnitude more efficient than brute-force particle-based 
algorithms il3| and this has allowed us to simulate gene 
expression on the relevant biological time scales of hours. 

Several publications [13, IH 111 111 EO, ISl El have 
discussed the effect of fluctuations in the binding of tran- 
scription factors to their site on the DNA (called opera- 
tor) on the noise in gene expression. Most of these models 
are relatively sim p le, i gno ring, for instance, production 
of mRNA 111 111, IT9II2I. Moreover all these studies, 
with the exception of [2^ [23 , ignore the role of the spa- 
tial fluctuations of the transcription factors. Our aim 
is to study gene expression in a biologically meaningful 
model. We have therefore constructed a rather detailed 
model, although we will also use minimal models that 
can be studied analytically, in order to interpret the sim- 
ulation results. The full model, which is described in the 
next section, contains the diffusive motion of repressor 
molecules, open complex formation, promoter clearance, 
transcription elongation and translation J24| . 



In section llVl we discuss the simulation results for both 
the noise in mRNA and protein level. The results reveal 
that for [R] < 50nM, the noise in the spatially-resolved 
model can be more than five times larger than the noise 
in the well-stirred model. We also show that a cell could 
minimize the effect of spatial fluctuations, either by tun- 
ing the open complex formation rate or by changing the 
number of repressors and their affinity for the binding 
site on the DNA. In section we elucidate the origin 
of the enhanced noise in the spatially resolved model. 
In the subsequent section, we show that in the model 
employed here the effect of spatial fluctuations can be 
quantitatively described by a well-stirred model in which 
the reaction rates for repressor binding and unbinding 
are appropriately renormalized; however, as we discuss in 
the last section, we expect that in a more refined model 
the effect of diffusion will be more complex, impeding 
such a simplified description. In section IVIII we discuss 
how the operator state fluctuations propagate through 
the different stages of gene expression using power spec- 
tra for the operator state, elongation complex, mRNA 
and protein. The results show that these power spectra 
are highly useful for unraveling the dynamics of gene ex- 



pression. We hope that this stimulates experimentalists 
to measure power spectra of not only mRNA and protein 
levels [23 , but also of the dynamics of transcription ini- 
tiation and elongation using e.g. magnetic tweezers '265. 
As we argue in the last section, such experiments should 
make it possible to determine the importance of spatial 
fluctuations for the noise in gene expression. 



II. MODEL 
A. Diffusive motion of repressors 

We explicitly simulate the diffusive motion of the re- 
pressor molecules in space. However, since the exper- 
iments of Riggs et al. [23 and the theoretical work of 
Berg, Winter, and Von Hippel |23|, it is well known that 
proteins could find their target sites via a combination 
of ID sliding along the DNA and 3D diffusion through 
the cytoplasm - "hopping" or "jumping" from one site 
on the DNA to another. This mechanism could speed 
up the search process and make it faster than the rate 
at which particles find their target by free 3D diffusion; 
this rate is given by fc = AnaD^ [R] , where a is the cross 
section, which is on the order of a protein diameter or 
DNA diameter, D3 is the diffusion constant of the pro- 
tein in the cytoplasm, and [R] is the concentration of the 
(repressor) protein. However, while it is clear that the 
mechanism of 3D diffusion and ID sliding could poten- 
tially speed up the search process, whether this mecha- 
nism in living cells indeed drastically reduces the search 
time is still under debate |29'|. In this context, it is in- 
structive to discuss the two main results of recent studies 
on this topic jM S IM HI IM & The first is that the 
mean search time r is given by |34| 
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where L is the total length of the DNA, A is the aver- 
age distance over which the protein slides along the DNA 
before it dissociates, Di is the diffusion constant for slid- 
ing, r is the typical mesh size in the nucleoid, and -D3 
is the diffusion constant in the cytoplasm. This formula 
has a clear interpretation [3J|: A^/Di is the sliding time, 
r^/Z?3 is the time spent on 3D diffusion, the sum of these 
terms is thus the time to perform one round of sliding and 
diffusion, and L/X is the total number of rounds needed 
to find the target. The other principal result is that the 
search time is minimized when the sliding distance A is 



(2) 



Under these conditions, a protein spends equal amounts 
of time on 3D diffusion and ID sliding (a protein is thus 
half the time bound to the DNA). Eq. |21is a useful re- 
sult, because it shows that the average sliding distance 
A depends upon the ratio of diffusion constants and on 




the typical mesh size in the nucleoid. If we now assume 
that Di and D^ are equal (which is not obvious given that 
proteins bind relatively strongly to DNA - Di could thus 
very well be much smaller than D3) and if we take the 
mesh size to be given by r ~ y/v/L |34| . where v « 1/im^ 
is the volume of an Escherichia coli cell and L w 10'^ /im, 
we find that A ~ lOnm (30 bp). This corresponds to 
the typical diameter of a protein or DNA double helix 
and is thus not very large. Interestingly, recent exper- 
iments seem to confirm this: experiments from Halford 
et al. on restriction enzymes (EcoRV and BbcCI) with a 
series of DNA substrates with two target sites and vary- 
ing lengths of DNA between the two sites, suggest that 
under the in vivo conditions, sliding is indeed limited to 
relatively short distances, i.e. to distances less than 50 
bp (« 16nm) [alll^- 

Now, it should be realized that on length scales be- 
yond the sliding length, the motion is essentially 3D dif- 
fusion: the sliding/hopping mechanism corresponds to 
3D diffusion with a jump distance given by the sliding 
distance 1^3 • Moreover, since the sliding distance is only 
on the order of a particle diameter, as discussed above, 
we have therefore decided to model the motion of the 
repressor molecules as 3D diffusion. But it should be re- 
membered that on length scales smaller than 10 — 30nm, 
this approach is not correct. We discuss the implications 
of this for our results in the Discussion section. 



B. Transcription and Translation 

Most repressors bind to a site that (partially) over- 
laps with the core promoter - the binding site of the 
RNA polymerase (RNAP). When a repressor molecule 
is bound to its operator site, it prevents RNAP from 
binding to the promoter, thereby switching off gene ex- 
pression. Only in the absence of a repressor on the op- 
erator site, can RNAP bind to the promoter and initiate 
transcription and translation, ultimately resulting in the 
production of a protein. We model this by the following 
reaction network: 
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Eqs. 13 and 2] describe the competition between the bind- 
ing of the repressor R and the RNAP molecules Rp to 
the promoter (O is the operator site). In our simula- 
tion we fix the binding site O in the center of a con- 
tainer with volume V — 1/im'^, comparable to the vol- 
ume of a single E. coli cell. We simulate both the op- 
erator site O and the repressor molecules as spherical 
particles with diameter a = lOnm. The operator site 
O is surrounded by Nn repressor molecules that move 
by free 3D diffusion (see previous section) with an ef- 
fective diffusion constant D = 1/im^s^^, as has been 
reported for proteins of a similar size |37| |. The intrin- 
sic forward rate km = 6 • 10^M~-^s~-^ for the repressor 
particles R at contact is estimated from the Maxwell- 
Boltzmann distribution hI^. The backward rate fcbR, de- 
pends on the interaction between the DNA binding site 
of the repressor and the operator site on the DNA and 
varies greatly between different operons, with stronger re- 
pressors having a lower /cbR- In our simulations, we vary 
/cbR between 1 — 0.01 s~^, as discussed in more detail 
below. The concentration of RNAP is much higher than 
that of the repressor [33l- Because of this we treat the 
RNAP as distributed homogeneously within the cell and 
we do not to take diffusion of RNAP into account explic- 
itly. Instead, RNAP associates with the promoter with 
a diffusion- limited rate k^p — 4:TTaD[Rp\. In our simula- 
tions, the concentration of free RNAP is [Rp] — 0.5/iAf 
|3q|. leading to a forward rate fcfRp = 38s^^. Finally, 
the backward rate /cbRp = 0.5 is determined such that 
Keci = 47rcrL>/fcbRp = 1-4 • lO^M"! H. 

Transcription initiation is described by Eqs. |S1 and 
El Before productive synthesis of RNA occurs, first the 
RNAP in the RNAP-promoter complex ORp unwinds ap- 
proximately one turn of the promoter DNA to form the 
open complex ORp*. The open complex formation rate 
koc has been measured to be on the order of 0.3 — 3s~^ 
[26j . We approximate open complex formation as an irre- 
versible reaction. Some experiments find this step to be 
weakly reversible J26i |. However, adding a backward reac- 
tion to the model did not change the dynamics of the sys- 
tem in a qualitative way, as long as the backward rate is 
smaller than fcoc j which is in agreement with experimen- 
tal results. After open complex formation, RNAP must 
first escape the promoter region before another RNAP or 
repressor can bind. Since elongation occurs at a rate of 
50 — 100 nucleotides per second and between 30 — 60 nu- 
cleotides must be cleared by RNAP before the promoter 
is accessible, a waiting time of tdcar = Is is required be- 
fore another binding can occur. Since promoter clearance 
consists of many individual elongation events that obey 
Poisson statistics individually, we model the step as one 
with a fixed time delay tcioar7 not as a Poisson process 
with rate 1/tcicar- 

Eqs. I7I11I describe the dynamics of mRNA and pro- 
tein numbers. After clearing the promoter region, RNAP 
starts elongation of the transcript T. As for clearance, 
the elongation step is modeled as a process with a fixed 
time delay ioion = 30s, corresponding to an elongation 



rate of 50 — 100 nucleotides per second and a 1500 bp 
gene. When a niRNA M is formed, it can degrade with 
a rate k^m- Here, the niRNA degradation rate is deter- 
mined by fixing the average mRNA concentration in the 
unrepressed state, as described below. Furthermore, a 
mRNA molecule can form a mRNA-ribosome complex 
Mribo and start translation. We assume that 6 = 5 
proteins are produced on average from a single mRNA 
molecule 3 , so that the start of translation occurs at a 
rate k^iho = b fcdm- After a fixed time delay itrans = 30s a 
protein P is produced. The mRNA is available for ribo- 
some binding immediately after the start of translation. 
Due to the delay in protein production, M can start to 
be degraded, while the mRNA-ribosome complex Mribo 
is still present; M thus represents the mRNA leader re- 
gion rather than the entire mRNA molecule. Finally, the 
protein P degrades at a rate fcdp, which is determined 
by the requirement that the average protein concentra- 
tion in the unrepressed state has a desired value, as we 
describe now. 

We vary the free parameters in the reaction network 
described in Eqs. 1311 II - Nr, fcbR, fcdm, fcdp - in the fol- 
lowing way: first, we choose the concentration of mRNA 
and protein in the absence of repressor molecules. In this 
case, tuning of the concentrations is most straightforward 
by adjustment of the mRNA and protein decay rates fcdm 
and fcdp. For the above reaction network one can show 
that the average mRNA number A^Af and protein number 
Np is given by 
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where Ki = kfRp/{kbpp + fcoc), ^2 = fcm/fcbR, K3 = 
fcoc^cioar, K4 = fcoc/fcdm and K5 = fcribo/fcdp are equi- 
librium constants, V is the volume of the cell and Np 
is the total number of repressors. The unrepressed state 
corresponds to Np = 0. In our simulations, we fix the 
mRNA and protein numbers in the unrepressed state at 
Nm = 50 and Np = 2-10^. The mRNA and protein de- 
cay rates then follow straightforwardly from Eqs. 1121 and 
nH the mRNA degradation rate is fcdm = 0.019s~i ^ 
and the protein degradation rate is fcdp = 2.4 x 10~'*s~"'^; 
the latter corresponds to protein degradation by dilution 
with a cell cycle time of around Ih. 

Next, we determine by what factor these concentra- 
tions should decrease in the repressed state. This can 
be done by changing the number of repressors Np and 
the repressor backward rate fcbR. We define the repres- 
sion level / as the transcription initiation rate in the ab- 
sence of repressors, divided by the initiation rate in the 
repressed state |41j. For a repression level /, the con- 
centration of mRNA and proteins in the repressed state 
is a fraction 1// of the concentration in the unrepressed 
state and it follows that 
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Thus, a fixed repression level / does not specify a unique 
combination of Np and fcbR: increasing the number of re- 
pressors twofold, while also increasing the repressor back- 
ward rate by the same factor, gives the same repression 
level. This means that the cell can control mRNA and 
protein levels in the repressed state either by having a 
large number of repressors that stay on the DNA for a 
short time or by having a small number of repressors, 
possibly even one, that stay on the DNA for a long time. 
Even though it is conceivable that the latter is preferable 
for economic reasons, there is no difference between the 
two extremes in terms of the average gene expression. In 
our simulations, we vary Np and fcbR, but use a fixed 
repression level / = 100. Consequently, in the repressed 
state, on average A^m = 0.5 and A'p = 200. 



III. SIMULATION TECHNIQUE 

We simulate the above reaction network using Green's 
Function Reaction Dynamics (GFRD) 0,0. As dis- 
cussed above, only the operator site O and the repressor 
particles R are simulated in space. All other reactions are 
assumed to occur homogeneously within the cell and are 
simulated according to the well-stirred model 42] or with 
fixed time delays for reaction steps involving elongation. 
A few modifications with respect to the algorithm de- 
scribed in |l4Lll5| are implemented to improve simulation 
speed. First, we neglect excluded volume interactions be- 
tween repressor particles mutually, as the concentration 
of repressor is very low. This means that the only po- 
tential reaction pairs we consider are operator-repressor 
pairs. Secondly, we use periodic boundary conditions in- 
stead of a reflecting boundary, which leads to a larger 
average time step. As the operator site O is both small 
compared to the volume of the cell and is far removed 
from the cell boundary, this has no effect on the dynam- 
ics of the system. Finally, as the repressor backward rate 
fcbR is rather small, the operator site can be occupied by 
a repressor for a time long compared to the average simu- 
lation time step. If the repressor is bound to the operator 
site longer than a time L^/6D, where L is the length of 
the sides of our container, the other repressor molecules 
diffuse on average from one side of the box to the other. 
Consequently, when the repressor eventually dissociates 
from the operator site, the other repressor molecules have 
lost all memory of their positions at the time of repressor 
binding. Here, when a repressor will dissociate after a 
time longer than L^ /6D, we do not propagate the other 
repressors with GFRD, but we only update the master 
equation and fixed delay reactions. We update the posi- 
tions of the free repressors at the moment that the oper- 
ator site becomes accessible again, by assigning each free 
repressor molecule a random position in the container; 
the dissociated repressor is put at contact with the op- 
erator site. We see no noticeable difference between this 
scheme and results obtained by the full GFRD algorithm 
described in Refs. [TUll. 
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FIG. 1; Dynamics of mRNA and protein numbers in the re- 
pressed state for different number of repressors Nr. The num- 
ber of mRNA and protein molecules is shown for simulations 
with GFRD (black line) and according to the master equation 
(gray line). In the GFRD simulation, diffusion of repressor 
particles is explicitly included, (a) and (b) Nr — 5. (c) and 
(d) Nr = 20. (e) and (f) Nr = 80. In general, there is a dra- 
matic difference in dynamics due to the spatial fluctuations 
of the repressor molecules. This difference becomes more pro- 
nounced as the number of repressors decreases. However, we 
find that in all cases (A^m) = 0.5 and (Np) = 200, on average. 



IV. SIMULATION RESULTS: DYNAMICS AND 
NOISE 



To study the effect of spatial fluctuations on the repres- 
sion of genes, we simulate the reaction network described 
in Eqs. l3llll both by GFRD, thus exphcitly taking into 
account the diffusive motion of the repressor particles, 
and according to the weU-stirred model, where the re- 
pressor particles are assumed to be homogeneously dis- 
tributed in space and the dynamics depends only on the 
concentration of repressor. In Fig. ^we show the behav- 
ior of mRNA and protein numbers for a system with open 
complex formation rate fcoc = 30s^^ and with varying 
numbers of repressors N^. We keep the repression factor 
fixed at / = 100 so that with increasing Nn the repres- 
sor backward rate fcbR is also increased, i.e. repressor 
particles are bound to the DNA for a shorter time. 

It is clear from Fig. ^ that there is a dramatic 
difference between the behavior of mRNA and protein 
numbers between the GFRD simulation and the well- 
stirred model. When spatial fluctuations of the repressor 
molecules are included, mRNA is no longer produced in 
a continuous fashion, but instead in sharp, discontinu- 
ous bursts during which the mRNA level can reach lev- 
els comparing to those of the unrepressed state. These 
bursts in mRNA production consequently lead to peaks 
in protein number. As the protein decay rate is much 
lower than that of mRNA, these peaks are followed by pe- 




FIG. 2: Noise in (a) mRNA number and (b) protein number 
as a function of the number of repressors A''^ and for constant 
repression factor / = 100. Data obtained by GFRD simula- 
tion is shown for fcoc = 0.3(o),3(n) and 30(*)s~^. Noise 
levels for the well-stirred model are shown as grey lines and 
those for the well-stirred model with reaction rates renormal- 
ized according to Eqs. 1161 and 1171 are shown as black lines, 
both for fcoc = 0.3 (solid lines), 3 (dashed lines) and 30 (dot- 
ted lines) s~^ . Only when the reaction rates are properly 
renormalized does the noise in the well-stirred model agree 
well with the noise in the GFRD simulations, which include 
the effect of diffusion. (Insets) Noise levels as a function of 
fcoc. Symbols indicate results for GFRD and lines are results 
for the chemical master equation with renormalized reaction 
rates. 



riods of exponential decay over the course of hours. Due 
to these fluctuations, protein numbers often reach levels 
of around 5 — 10% of the protein levels in the unrepressed 
state. In contrast, in the absence of repressor diffusion, 
the fluctuations around the average protein number are 
much lower. For both cases, however, the average behav- 
ior is identical: even though the dynamics is very differ- 
ent, we always find that on average (iVmRNA) = 0.5 and 
(Np) = 200. Also, in all cases the fluctuations in mRNA 
number are larger than those in protein number. This 
means that the translation step functions as a low-pass 
fllter to the repressor signal. 

When we increase the number of repressors Np and 
change fcbR in such a way that the repression level / re- 
mains constant, we find that both for GFRD and the 
well-stirred model the fluctuations in mRNA and pro- 
tein number decrease. In the absence of spatial fluctu- 
ations this effect is minor, but for GFRD this decrease 
is sharp: for large number of repressors, the burst in 
mRNA become both weaker and more frequent. This in 
turn leads to smaller peaks and shorter periods of expo- 
nential decay in protein numbers. In fact, as Np is in- 
creased both approaches converge to the same behavior. 
At around Np « 100, the dynamics of the protein num- 
ber is similar for the well-stirred model and the spatially 
resolved model. The same happens for mRNA number 
when Np ^ 500. In Fig. |21 we quantify the noise in 
mRNA and protein number, deflned as standard devia- 
tion divided by the mean, while we change the number 
of repressors Np. As we keep the amount of repression 
flxed at / — 100, we simultaneously vary the backward 
rate fcbR according to Eq. ^1 When all parameters are 



the same, the noise for the GFRD simulation, including 
the diffusive motion of the repressors, is always larger 
than the noise for the well-stirred model, where the diffu- 
sive motion is ignored. In both cases, the noise decreases 
when the number of repressors is increased and the re- 
pressor backward rate becomes larger. This is consistent 
with the niRNA and protein tracks shown in Fig. ^ We 
also investigated the effect of changing the open complex 
formation rate koc- In nature, this rate can be tuned 
by changing the base pair composition of the promoter 
region on the DNA. When we change koci we change 
the mRNA decay rate fcdm so that the average mRNA 
and protein concentrations remain unchanged (see sec- 
tion lIIB|l . We find that when hoc is lowered, the fluctu- 
ations in mRNA and protein levels are sharply reduced. 
When koc is much larger than the RNAP backward rate 
^bRp = 0.5s~^, almost every RNAP binding to the pro- 
moter DNA will result in transcription of a mRNA. For 
/cqc smaller than fcbRpi RNAP binding will lead to tran- 
scription only infrequently. As a consequence, the oper- 
ator filters out part of the fluctuations in RNAP bind- 
ing due to the diffusive motion of the repressor particles, 
leading to the decrease in noise observed in Fig. |21 This 
shows that the open complex formation rate plays a con- 
siderable role in controlling noise in gene expression. 



V. SIMULATIONS RESULTS: OPERATOR 
BINDING 

To understand how the diffusive motion of repressor 
molecules leads to increased fluctuations in mRNA and 
protein numbers, it is useful to look in some detail at the 
dynamics of repressor-DNA binding. In figure I^V, we 
show the OR bias for both GFRD and the well-stirred 
model. The OR bias is a moving time average over OR{t) 
with a 50s time window and should be interpreted as the 
fraction of time the operator site was bound by repres- 
sor particles over the last 50 seconds. The results we 
show here are for Nr = 5 repressors and a repression 
factor / = 2. At this repression factor, fcbR is such that 
the repressor molecules are bound to the operator only 
fifty percent of the time, making it easier to visualize the 
operator dynamics than in the case of / = 100 as used 
above. The OR bias for the wcU-stirred model fiuctuates 
around the average value {OR) — 0.5, indicating that on 
the timescale of 50s several binding and unbinding events 
occur, in agreement with fcbR, = 1.26s~^ for / = 2. On 
the other hand, when including spatial fiuctuations, the 
OR bias switches between periods in which repressors are 
bound to the DNA continuously and periods in which 
the repressors are virtually absent, both on timescales 
much longer than the 50s time window. How is it pos- 
sible that repressors are bound to the operator site for 
times much longer than the timescale set by the dissoci- 
ation rate from the DNA? The answer to that question 
can be found in Figs. |2J3 and C, where a time trace 
is shown of the operator occupancy by the repressor for 
both GFRD and the well-stirred model. The time trace 




FIG. 3: Dynamics of repressor binding for a repression factor 
of / = 2 and Nr = 5. (a) The Oi?-bias for GFRD (black 
line) and the well-stirred model (gray line). The O-R-bias 
is defined as the fraction of time a repressor is bound to the 
operator site in the last 50 seconds. When the diffusive motion 
of repressor molecules is included (black line), the O-R-bias 
switches between periods where repressors are continuously 
bound to or absent from the DNA for long times, (b) and (c) 
Time trace of the occupancy of the operator site by repressor 
molecules. When OR = 1 a repressor is bound to the operator 
site and OR = indicates either a free operator site or one 
with RNAP bound. For the GFRD simulations, an initial 
binding is followed by several rapid rebindings, whereas for 
the well-stirred model binding and rebinding is much more 
unstructured. Note that here, for reasons of clarity, f — 2 
instead oi f — 100 as used in the text and Figs. 0and|5| 



for the simulation of the well-stirred model in Fig. |3P 
shows a familiar picture: binding and dissociation of the 
repressor from the operator occurs irregularly, the time 
between events given by Poisson distributions. The time 
trace for GFRD in Fig. 03 looks rather different. Here, 
in general a dissociation event is followed by a rebind- 
ing very rapidly. Only occasionally does a dissociation 
result in the operator being unbound by repressors for 
a longer time. When this happens, repressors stay away 
from the operator for a time much longer than the typical 
time separating binding events in Fig. |2p. These series 
of rapid rebindings followed by periods of prolonged ab- 
sence from the operator result in the aberrant OR bias 
shown in Fig. |2K- 

The occurrence of rapid rebindings is intimately re- 
lated to the nature of diffusion. When diffusion and 
the positions of the reactants are ignored all dynamics is 
based only on the average concentration of the reactants. 
As a consequence, when in this approach a repressor dis- 
sociates from the operator site, the probability of rebind- 
ing depends only on the concentration of repressor in the 
cell. On the level of actual positions of the reactants, 
this amounts to placing the repressor at a random posi- 
tion in the container. The situation is very different for 
the GFRD approach, where the positions of the reactants 
are taken into account. After a dissociation from the 
operator site, the repressor particle is placed at contact 
with the operator site. Because of the close proximity of 
the repressor to its binding site, it has a high probabil- 
ity of rapidly rebinding to, and only a small probability 



of diffusing away from, the binding site. At the same 
time, when the repressor eventually diffuses away from 
the operator site, the probability that the same, or more 
likely, another repressor diffuses to and binds the opera- 
tor site is much smaller than the probability of binding 
in the well-stirred model, as will be shown quantitatively 
in Sec. IVII This results in the behavior observed in Fig. 

It can now be understood that the bursts in mRNA 
production correspond to the prolonged absence of re- 
pressor from the operator site compared to the well- 
stirred model. Especially for low repressor concentra- 
tions, these periods of absence can be long enough that 
the concentration of mRNA reaches values comparable to 
those in the unrepressed state for brief periods of time. 
When a repressor binds to the operator site, due to the 
rapid rebindings it will remain bound effectively for a 
time much longer than the mRNA lifetime, leading to 
long periods where mRNA is absent in the cell. This 
shows that under these conditions spatial fluctuations 
and not stochastic chemical kinetics are the dominant 
contribution to the noise in mRNA and protein numbers 
in the repressed state. 



VI. TWO-STEP KINETIC SCHEME 

In this section we investigate to what extent the effect 
of diffusion on the repressor dynamics can be modeled by 
the two-step kinetic scheme [43, Hj '■ 



+ R^O---R^OR. 



(15) 



The first step in Eq. El describes the diffusion of re- 
pressor to the operator site resulting in the encounter 
complex O ■ ■ ■ R, with the rates fc_|_ and fc_ depending 
on the diffusion coefficient D and the size of the parti- 
cles. The next step describes the subsequent binding of 
repressor to the DNA. In this case the rates are related 
to the microscopic rates defined in Eq. |31 When the 
encounter complex is assumed to be in steady state, the 
two-step kinetic scheme can be mapped onto the reac- 
tion described in Eq. O but with effective rate constants 
k'fj^ = k+ka/{k- + ka) and k'^^ = k_kd/{k_ + ka) O. 
The two-step kinetic scheme should yield the same aver- 
age concentrations as the scheme in Eq. El so that the 
equilibrium constant K = ka/kd = ^fi^/^hi^ = kfn/khi^, 
where fcfR and fcbR are the reaction rates defined in Eq. 

m 

It is possible to express the effective rate constants 
/c^ and fc^jfj^ in terms of the microscopic rate constants 
kfR and khR- For the setup used here, where a single 
operator O is surrounded by a homogeneous distribution 
of repressor R, the rate k+ follows from the solution of the 
steady state diffusion equation with a reactive boundary 
condition with rate k = ka a,t contact p4l |45| and is 
given by the diffusion-limited reaction rate ko = ^ttctD. 
The rates fc_ and ka depend on the exact definition of 



the encounter complex O ■ ■ ■ R. It is natural to identify 
the rate kd with the intrinsic dissociation rate kbR, thus 
kd — kbR. From these expressions for fc+ and kd and the 
requirement that the equilibrium constant should remain 
unchanged, one finds that ka/k- — kfR/kr). Using this 



result one obtains k'^^ 



kokfR/iko + kfR) and fc^ 



bR 



kokbR/iko + kfR). 

These renormalized rate constants have a clear inter- 
pretation. For the effective forward rate it follows, for 
instance, that: \/k'r^ = l/ku + l/kfR-. that is, on av- 
erage, the time required for repressor binding is given 
by the time needed to diffuse towards the operator plus 
the time for a reaction to occur when the repressor is in 
contact with the operator site jijl- The effective back- 
ward rate has a similar interpretation. The probabil- 
ity that after dissociation the repressor diffuses away 
from the operator site and never returns is given by 
Siri{t -^ oo|(t), where S'iri(i, ''o) is the irreversible sur- 
vival probability for two reacting particles |4y]. Using 
that Siri{t -^ co\a) — k]j/{kfR + k^), the expression for 
k'fjj^ can be written as A;^^ = kbRSirr{t — > oo|cr): that is, 
the effective dissociation rate is the microscopic dissoci- 
ation rate multiplied by the probability that after disso- 
ciation the repressor escapes from the operator site |44J| . 

For diffusion limited reactions, such as the reaction 
considered here, we have that kfR ^ fc/^. Now, the 
renormalized rate constants reduce to: 



^fR — ^D, 

Kr = kokbR/kfR. 



(16) 
(17) 



In Fig. 121 we compare the noise profiles for the GFRD 
algorithm with those obtained by a simulation of the well- 
stirred model, where instead of the microscopic rates kfR 
and kfB we use the renormalized rates from Eqs. 1161 
and 1171 Surprisingly, we find complete agreement. One 
of the main reasons why this is unexpected, is that for 
the master equation the time between events is Poisson- 
distributed, whereas after a dissociation the time to the 
next rebinding is distributed according to a power-law 
distribution when diffusion is taken into account 4(|. 
The reason that this power-law behavior of rebinding 
times is not of influence on the noise profile, is that the 
time scale of rapid rebinding is much smaller than any 
of the other relevant time scales in the network. Specif- 
ically, rebinding times are so short that the probability 
that a RNAP will bind before a rebinding is negligible. 
As a consequence, the transcription network is not at 
all infiuenced by the brief period the operator site is ac- 
cessible before rebinding: for the transcription machinery 
the series of consecutive rebindings, albeit distributed al- 
gebraically in time individually, is perceived as a single 
event. And on much longer time scales, when a repres- 
sor diffuses in from the bulk towards the operator site, 
the distribution of arrival times is expected to be Pois- 
sonian, because on these time scales the repressors are 
distributed homogeneously in the bulk. 

It is possible to reinterpret the effective rate constants 
in Eq. 1161 and 1171 in the language of rapid rebindings. 
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The probability p that a rebind will occur after a disso- 
ciation ft'om the DNA is given hy p — 1 — Soo, where 
St — Sii-i{t,rQ = cr). The probability that n consecutive 
rebindings occur before the repressor diffuses away from 
the operator site is then given by p„ = (1 — SocY^Sqc. 
From this follows that the average number of rebind- 
ings is Nrb = (l - Soo)/Soa- Using again that S'oo == 
ko/ikfR -t- fcu), we find that Nub = kfn/kn- Combin- 
ing this with Eqs. 1161 and [T7I we get: 



k'fR 
Hr 



kfR/NRB, 

kbR/NRB- 



(18) 
(19) 



In words, after an initial binding the repressor spends 
Nrb times longer on the DNA than expected on the ba- 
sis of the microscopic backward rate, as it rebinds on av- 
erage Nrb times. Because the average occupancy should 
not change, the forward rate should be renormalized in 
the same way. In conclusion, in this model the effects 
of diffusion can be properly described by a well-stirred 
model when the reaction rates are renormalized by the 
average number of rebindings. 



VII. POWER SPECTRA 

In this section, we study how the noise due to the 
stochastic dynamics of the repressor molecules propa- 
gates through the different steps of gene expression for 
both the spatially resolved model and the well-stirred 
model. This analysis will also provide further insight 
into why the well-stirred model with renormalized rate 
constants for the (un) binding of the repressor molecules 
works so well. 

In biochemical networks, the noise in the output signal 
depends upon the noise in the biochemical reactions that 
constitute the network, the so-called intrinsic noise, and 
on the noise in the input signal, called extrinsic noise 
II m H El E3- In our case, the output signal is 
the protein concentration, while the input signal is pro- 
vided by the repressor concentration. The intrinsic noise 
arises from the biochemical reactions that constitute the 
transcription and translation steps. Moreover, we con- 
sider the noise in the protein concentration that is due 
to the (un)binding of the RNAP to (from) the DNA to 
be part of the intrinsic noise. The extrinsic noise is pro- 
vided by the fluctuations in the binding of the repressor 
to the operator, i.e. in the state OR. Since the total 
repressor concentration, [Rt\ ~ [R\ + [Oi?], is constant, 
the extrinsic noise is also given by the fluctuations in the 
concentration of unbound repressor. 

The noise properties of biochemical networks are most 
clearly elucidated via the power spectra of the time traces 
of the copy numbers of the components. Recently, we 
have shown that if the fluctuations in the input signal 
are uncorrelated with the noise in the biochemical reac- 
tions that constitute the processing network the power 
spectrum of the output signal is given by [50j 



S'p(cj) = S'in(t^) + g{uj)Scy,{uj). 



(20) 



Here, 5'p(w) is the power spectrum of the output signal, 
the protein concentration. The spectrum Sin{(^) denotes 
the intrinsic noise of the processing network; it is defined 
as the noise in the output signal in the absence of noise 
in the input signal. Here, the intrinsic noise is due to 
the biochemical reactions of transcription and transla- 
tion. The spectrum Scx{^) is the power spectrum of the 
input signal, which, in this case, is given by the noise in 
the concentration of unbound repressor: S'cx('^) = Sb.{^)] 
because the total repressor concentration is constant this 
power spectrum is also directly related to that of the 
repressor-bound state of the operator, S'or(i^). The func- 
tion g{uj) is a transfer function, which indicates how fiuc- 
tuations in the input signal are transmitted towards the 
output signal. If the extrinsic noise is uncorrelated with 
the intrinsic noise, then g{Lo) is an intrinsic quantity that 
only depends upon properties of the processin g ne twork, 
and not upon properties of the incoming signal 5Ci| . How- 
ever, for the network studied here, the noise in the input 
signal is not uncorrelated with the intrinsic noise ||5C1| . 
As we have shown recently, this means that Eq. [201 is 
not strictly valid |50|; the extrinsic contribution to the 
power spectrum of the output signal can no longer be 
factorized into a function that only depends upon intrin- 
sic properties of the network, g{Lo), and one that only 
depends upon the input signal, Sci^{lj). This relation is 
nevertheless highly instructive. Indeed, Eg. 1201 could be 
interpreted as a heuristic definition of the transfer func- 
tion g(u}). 

The diffusive motion of the repressor molecules impede 
an analytical evaluation of the power spectrum for the 
extrinsic noise. Moreover, while power spectra can be 
calculated analytically for linear reaction networks |5ll| . 
the delays in transcription resulting from promoter clear- 
ance and elongation, preclude the derivation of an ana- 
lytical expression for the power spectrum of the intrin- 
sic noise. We have therefore obtained the power spectra 
S-piuj), Sc-k{uj), and 5'in(w), directly from the time traces 
of the copy numbers. The power spectrum of a compo- 
nent X is given by S'x(w) — {\X{lo)\'^), where X{lo) is 
the Fourier Transform of the concentration X{t) of com- 
ponent X. Conventional FFT algorithms are not conve- 
nient, because our signals vary over a wide range of time 
scales. We therefore adopted a novel and efficient proce- 
dure, which is described in the appendix. This procedure 
should prove useful for computing the power spectra of 
time traces of copy numbers of species in biochemical net- 
works, as obtained by Kinetic Monte Carlo simulations. 

As indicated above, the intrinsic noise, Si-a_{Ljj), is de- 
fined as the noise in the output signal in the absence of 
fiuctuations in the input signal. In order to determine 
the intrinsic contribution to the noise in the protein con- 
centration, we discarded the (un)binding reaction of the 
repressor to the DNA (Eq. O, while rescaling the rate 
fcbRp for the dissociation reaction of the RNAP from the 
DNA (Eq. ^ in such a way that the average concentra- 
tion of the protein P remains unchanged. This eliminates 
the extrinsic noise arising from the repressor dynamics. 



thereby allowing us to obtain the intrinsic noise of the 
reactions in Eqs. 1411 II The rescaled backward rate fc^j^p 
is given by 



'^bRp 



fcbRp(l + K2NR/V) + kocK2/V 



(21) 



where K2 = fcfR/fcbR. 

For the interpretation of the power spectra of the 
niRNA and protein concentration, as discussed below, 
it is instructive to recall the power spectrum of a linear 
birth-and-death process, 



4aA 



(22) 



with rate constants k and fi. For the interpretation of 
the spectra of repressor binding to the DNA, it is useful 
to recall the spectrum of a two-state model. 



0^0*, 

fc2 



(23) 



with rate constants fci and ^2- For both models, the 
power spectrum is a Lorentzian function of the form: 



S{u;) 



20-2^ 

fjp + 07" 



(24) 



For the birth-and-death process, the variance in the con- 
centration of ^, cr^, is k/fi, while for the two state system, 
the variance a'^ in the occupancy n is n(l — n); the decay 
rate in the two-state model is /i = fci + /c2. The corner 
frequency fj, yields the time scale on which fluctuations 
relax back to steady state. We also note that the noise 
strength a^ is given by the integral of the power spectrum 
S{lu): a"^ = l/(27r)/^ dwS'(aj). The noise strength is 
thus dominated by those frequencies at which the power 
spectrum is largest. 

In the next subsection, we discuss the effect of spa- 
tial fluctuations on the noise in gene expression and ex- 
plain why a well-stirred model with renormalized rate 
constants for repressor (un)binding can capture its ef- 
fect. In the subsequent section, we discuss how the noise 
is propagated through the different stages of gene expres- 
sion. 



A. Spatial Fluctuations 

In Fig. ^ we show the power spectra for the input and 
output signals, for both the spatially resolved model and 
the well-stirred model with renormalized rate constants 
for repressor (un)binding (see previous section). We re- 
call that the output signal is the protein concentration, 
while the input signal is the concentration of unbound re- 
pressor (the extrinsic noise). Fig. ^ also shows the power 
spectrum of the intrinsic noise. This is the noise in the 
protein concentration (the output signal) , when the noise 
in the input signal resulting from the repressor dynam- 
ics has been eliminated by the procedure outlined above. 
The power spectra have been obtained in a parameter 
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FIG. 4: Power spectra of the repressor and protein concentra- 
tions obtained for / = 100, fcoc ~ 30s~^, Nr — 5. Data are 
shown both for the renormalized well-stirred model (RWS) 
with reaction rates renormalized according to Eqs. I16I17I 
and for GFRD, taking into account the spatial fluctuations 
of the repressor molecules explicitly. Also shown is the power 
spectrum of the intrinsic noise, which is the power spectrum 
of the protein concentration in the absence of fluctuations 
in the repressor concentration (extrinsic noise). For large uj, 
the repressor spectrum (extrinsic noise) differs between the 
well-stirred and the spatially resolved model. However, this 
difference does not appear in the power spectra of the protein 
concentration. The inset shows the frequency-dependent gain 
g{Lo) (see Eg.l^. 



regime where the diffusing repressors have a large effect 
on the noise: koc = 30s~^, Nji = 5 (see Fig.lSJ. 

Fig. ^ shows that the power spectrum of the protein 
concentration in the spatially resolved model is identical 
to that in the well-stirred model for the entire range of 
frequencies observed. This confirms the observation in 
Section IVII that the effect of the spatial fluctuations of 
the repressor molecules on the noise in the protein con- 
centration can by described by a well-stirred model in 
which the reaction rates for repressor (un)binding to the 
DNA are properly renormalized. 

Fig. ^ also elucidates the reason why a well-stirred 
model with properly renormalized rate constants for re- 
pressor (un)binding can successfully describe the effect 
of the diffusive motion of the repressor molecules on the 
noise in gene expression. It is seen that the repressor 
spectrum for the renormalized well-stirred model is ac- 
curately described by a Lorentzian function with a cor- 
ner frequency /i = 0.02s~^, as expected for the dynam- 
ics of repressor (un)binding dynamics (see next section). 
The repressor spectrum of the spatially resolved model 
fully overlaps with that of the well-stirred model up to 
a frequency of cj « lO^s^^, but for higher frequencies 



it shows a clear deviation from the 



behavior. This 



deviation is caused by the diffusive motion of the repres- 
sor molecules. Indeed, the deviation occurs at frequen- 
cies comparable to the inverse of the typical time scale 
for rapid rebindings {^ fis). However, this difference be- 
tween the spectrum of the repressor dynamics in the spa- 
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creases as the relaxation rate fj, = ki + k2 decreases as 
a result of the slower binding and unbinding of repressor 
(see Eq. 124(1 . The higher power in the repressor spec- 
trum at low frequencies for the spatially resolved model 
and for the well-stirred model with the renormalized rate 
constants, as compared to that for the well-stirred model 
with the unrenormalized rate constants, is not filtered 
by the processing network of transcription and transla- 
tion and thus manifests itself in the power spectrum of 
the protein concentration. Spatial fluctuations of gene 
regulatory proteins thus increase the noise in gene ex- 
pression by increasing the power of the input signal at 
low frequencies. 



FIG. 5: The power spectra for the well-stirred model 
with unrenormalized rate constants (WS) and for the well- 
stirred model with renormalized rate constants for repressor 
(un)binding to (from) the DNA (RWS). The intrinsic noise 
of gene expression is the same for both models. The ex- 
trinsic noise, arising from the repressor dynamics, is, how- 
ever, markedly different. The repressor spectrum for the well- 
stirred model with renormalized rate constants has lower cor- 
ner frequency, but, more importantly, also a higher power at 
low frequencies. The increased power at low frequencies is not 
filtered by the processing network and increases the noise in 
gene expression. 



tially resolved model and that in the well-stirred model 
does not manifest itself in the spectra for the protein con- 
centrations of the two respective models, for two reasons: 
1) the difference only occurs at high frequencies, i.e. in a 
frequency regime where the fluctuations only marginally 
contribute to the noise strength (the difference in area 
under the curves of the repressor power spectra for the 
two models is less than 5%); 2) the repressor fluctuations 
in this frequency range are filtered out by the processing 
network of transcription and translation; as a result of 
this, the effect of the small difference in area under the 
curves of the repressor power spectra for the two models 
is reduced even further. The filtering properties of the 
processing network are illustrated in the inset of Fig. 01 
which shows the transfer function g{uj) as obtained from 
g{uj) = (S'p(w) - S'in(t^))/S'ex(w) (sce Eq. mj. Clearly, 
the transfer function rapidly decreases as the frequency 
increases. This shows that the processing network of 
transcription and translation acts as a low-pass filter, re- 
jecting the high frequency noise in the repressor dynam- 
ics that originates from the rapid rebindings. The only 
effect of the repressor rebindings on the noise in gene ex- 
pression is thus that it lowers the effective dissociation 
rate (and association rate), as explained in the previous 
section. As compared to the well-stirred model with the 
unrenormalized rate constants for repressor (un)binding, 
this decreases the corner frequency /j, in the repressor 
power spectrum (see Fig. EJ, but increases the power 
at low frequencies - recall that for a two-state model, 
which relaxes mono-exponentially, the power spectrum 
at zero frequency is S{uj = 0) = 2cr^//i, which thus in- 



B. Noise propagation 

In Fig. we show how fluctuations in the input sig- 
nal arising from the dynamics of repressor binding and 
unbinding, are propagated through the different stages 
of gene expression. In Fig. El^a) we illustrate how the 
noise in the repressor concentration (the extrinsic noise) 
is transferred to the level of transcription. The figure 
shows for both the spatially resolved model and for the 
well-stirred model with renormalized rate constants for 
repressor (un) binding, the power spectrum of the repres- 
sor concentration and the spectrum of the concentration 
of the elongation complex, defined as [ORp*] + [T]. It is 
clear from Fig. IHIa) that already at the level of the elon- 
gation complex, the high-frequency noise due to the rapid 
rebindings is filtered. Transcription can thus already be 
described by a well-stirred model with properly renormal- 
ized rate constants for repressor (un) binding to (from) 
the DNA. The power spectrum of the elongation complex 
exhibits two corner frequencies, one around u;+ ~ 40s~^ 
and another one at w_ « 0.02s^^. These two corner 
frequencies arise from the competition between repressor 
and RNAP for binding to the promoter. To elucidate 
this, we have plotted in the inset the power spectrum for 
RNAP bound to the promoter, thus the power spectrum 
for [ORp] + [ORp*]. It is seen that this power spectrum 
has the same two corner frequencies as that of the elonga- 
tion complex, showing that their dynamics is dominated 
by the same processes - repressor binding and RNAP 
binding to the promoter. These two corner frequencies 
can be estimated analytically by considering the reac- 
tions in Eqs. I3I6I as a three-state system, in which re- 
pressor and RNAP compete for binding to the promoter: 



^2 ks 

OR^Ot± ORp'. 



(25) 



Here, ORp' — ORp + ORp*, where ORp denotes the 
RNAP bound to the promoter in the closed complex 
and ORp* denotes RNAP bound to the promoter in the 
open complex. The rate constant fci denotes the rate 
at which a repressor binds to the promoter; it is given 



by fci = k'^ijlRr] 



where k'^j^ is the renormalized asso- 
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FIG. 6: Comparison of the power spectra at different stages 
of gene expression, (a) Power spectrum for repressor con- 
centration and for the elongation complex ORp* + T, both 
for the well-stirred model with renormalized rate constants 
(RWS) and for GFRD. Repressor power spectra show a differ- 
ence between the spatially resolved model and the well-stirred 
model at high frequencies, due to the diffusion of the repressor 
molecules. The power spectra for the elongation complexes 
coincide for the well-stirred and the spatially resolved model. 
The power spectrum of the elongation complex shows a se- 
ries of peaks and valleys due to the presence of fixed delays 
in the dynamics of the elongation complex, (inset) Power 
spectrum of RNAP dynamics {ORp + ORp). Shown are the 
power spectra in the presence and absence of fixed delays in 
the RNAP dynamics. Due to the competition between RNAP 
and rerepssor for binding to the promoter, the power spec- 
trum is described by a sum of two Lorentzians. (b) Power 
spectra of the elongation complex and mRNA. Peaks due to 
the delays in RNAP dynamics are still present in the mRNA 
dynamics. For high frequencies, the mRNA dynamics is well 
described by a linear birth-and-death process, (c) Spectra of 
mRNA and protein. The slow protein dynamics filters out all 
the peaks resulting from the delays in the RNAP dynamics. 
The only difference between the full spectrum of the output 
signal and that of the intrinsic noise is an increased noise at 
low frequencies, due to the repressor dynamics. 



elation rate (see Eg. I16fl . The rate constant ^2 denotes 
the renormalized rate for repressor unbinding, ^2 = fc(,jj, 
(see Eq. I17|l ; k^ = kfnp denotes the rate at which RNAP 
binds to the promoter. The rate constant k^ is the rate 
at which the RNAP leaves the promoter. Since the pro- 
moter can become accessible for the binding of another 
RNAP or repressor by either the dissociation of RNAP 
from the closed complex or by forming the open com- 
plex and then clearing the promoter, this rate is given 
by ki = kbRp + {kp^ + icicar)^^- If promoter clearance 
would be neglected, then, indeed, k^ = kbup + koc- 

The power spectrum of the RNAP dynamics in Eq. 1251 
can be calculated analytically and is given by a sum of 



two Lorentzians: 

SoRp' (w) 



Auj- 



Buj, 



■"+ 



-uj^ 



(26) 



where A and B are coefficients. The corner frequencies 
LO- and 0;+ are given by uj± = (fc ± \Jk'^ — Ah)/2, where 
k — ^^ki and h — kik^ + k2{k^ + k^). The dynam- 
ics of repressor binding and unbinding is much slower 
than that of RNAP binding and unbinding, meaning that 
fci, fc2 <C fcs, A:4. This allows us to approximate the corner 
frequencies as a;+ = k^ + k^ and w_ = fc2 + fcifc4/(fc3 + fc4)- 
This yields the following expressions for the corner fre- 
quencies: 



^fRp 



+ kfjRp + [kpfj + ^clear) 



k'bR + k'[RT][0]'. 



(27) 
(28) 



Here, [O]' = ki/{k^ + k^) is the conditional probability 
that the promoter is not occupied by the RNAP, given 
that it is not occupied by repressor; it is given by the 
occupancy of the promoter by RNAP in the absence of 
any repressor molecules in the system. We can now see 
that the highest corner frequency, W-(_, describes the fast 
dynamics of RNAP binding to, and clearing from, the 
promoter and that the other corner frequency, lo^, rep- 
resents the slow dynamics of repressor (un)binding to the 
DNA in the presence of the fast RNAP bindings to the 
promoter; the lower corner frequency, W-, is also the cor- 
ner frequency in the repressor spectrum of the renormal- 
ized well-stirred model (see Figs. 0] and Isjl. In Fig. |B|[a) 
we plot the power spectrum SoRp' (w) as predicted by the 
three-state model (Eq. |5S1 with fitted coefficients A and 
B) on top of the power spectrum obtained from the sim- 
ulations and find excellent agreement. We also show the 
power spectra when we neglect the delay due to promoter 
clearance. As expected, in the absence of the delay due to 
promoter clearance, the lower corner frequency, lu-, and, 
to a smaller extent, the higher corner frequency, w+, are 
shifted to higher frequencies. 

The power spectrum of the elongation complex in Fig. 
EJa) contains information that is not easily observed in 
the time domain and could as a result be helpful in 
the interpretation of the results. It is seen that there 
are two series of peaks. Those are associated with the 
two processes with fixed time delays. The first pro- 
cess is the promoter clearance, which takes a fixed time 
tcicar- Indeed, the ffist peak in the corresponding series 
of peaks in the power spectrum of the elongation com- 
plex, is at a; ~ 27r/(icicar) = 6.3s~^; the other peaks 
in the series are the higher harmonics that naturally 
arise for processes with fixed time delays. The second 
process is the transcript elongation process. After the 
elongation complex has been formed, it takes a fixed 
time icicar + ioion beforc the full transcript is formed and 
the RNAP dissociates from the DNA; the first valley of 
the corresponding series of peaks/valleys is, indeed, at 
uj w 27r/(tcicar + ^cion) = 0.2s^^. While the frequency 
27r/icicar ylclds, to a good approximation, the rate at 
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which the elongation complex signal increases, the fre- 
quency 27r/(iciear + ieion) corresponds to the frequency at 
which the elongation complex signal decreases; this ex- 
plains why the shapes of the respective series of peaks 
and valleys are reciprocal. Lastly, the reason that both 
peaks and valleys are broadened is that the delay in the 
formation of the elongation complex is not fully determin- 
istic: the duration of the delay is not only determined by 
the promoter clearance time, which, indeed, is fixed, but 
also by the time it takes for another RNAP to bind the 
DNA and then form the open complex - in the absence of 
repressor, the average frequency at which an elongation 
complex is formed is given by 27r/(fcj^ -I- k^^ + icicar) 
(see also Eqs.0]-|5l. Both RNAP binding and open com- 
plex formation are modeled as Poisson processes, and this 
leads to a distribution of delay times for the formation of 
the elongation complex. 

In Fig. Elb) and (c) , we examine how the noise in the 
dynamics of the elongation complex propagates to the 
level of niRNA and protein dynamics. In Fig. EJb), we 
compare the full power spectrum of the mRNA concen- 
tration with that of the elongation complex - the input 
signal (extrinsic noise) for the mRNA signal - and that of 
the intrinsic noise of the mRNA signal; to compute the 
intrinsic noise, we have modeled the mRNA dynamics 
as a birth-and-death process (see Eq. I22|) with a pro- 
duction rate as given by the average production rate for 
the full system in Eqs. I3I11I As expected, for higher 
frequencies {u; > O.ls^^), the full spectrum of mRNA 
overlaps almost fully with that of the intrinsic noise, al- 
though some traces of the input signal (the elongation 
complex) are still apparent in this high frequency regime; 
these are the peaks at iv ~ 6.3s~^ corresponding to pro- 
moter clearance. At lower frequencies (w < 0.1s~^), the 
noise in the mRNA signal is dominated by the extrinsic 
noise, which is the noise in the elongation complex (the 
input signal). Indeed, both the spectrum of the elonga- 
tion complex and that of mRNA have a corner frequency 
at uj-, which, as discussed above, arises from the slow 
repressor (un)binding to the DNA in the presence of the 
fast DNA-(un)binding kinetics of RNAP. 

Fig.Efc) shows how the noise in the mRNA concentra- 
tion is propagated to that in the protein concentration. 
Again, at higher frequencies, the spectrum of the protein 
concentration coincides with that of the intrinsic noise of 
protein synthesis, which, as above for mRNA, has been 
computed by modeling protein production as a birth-and- 
death process; note also that the remnants of operator 
clearance (the peaks in the spectrum at w w 6.3s^^) have 
been filtered by the slow protein dynamics. Only for fre- 
quencies smaller than uj ~ 0.1s~^, does the extrinsic noise 
- the noise in the mRNA concentration - strongly con- 
tribute to the noise in the protein concentration. A care- 
ful inspection of the protein spectrum shows that it has 
a "corner" at aj_ , which arises from the repressor DNA- 
(un)binding dynamics (the extrinsic noise), and one, al- 
beit much less visible, at w w fcdp = 2 x 10^'^s^^, which 
is due to the intrinsic dynamics of protein degradation. 



VIII. DISCUSSION AND OUTLOOK 

Our analysis reveals that at high frequencies both 
mRNA and protein synthesis are well described by a lin- 
ear birth-and-death model. In this frequency regime, the 
effect of spatial fluctuations, originating from the rapid 
repressor rebindings, is completely filtered by the slow 
dynamics of transcription and translation. These rebind- 
ings do, however, decrease the effective rate at which the 
repressor molecules associate with, and dissociate from, 
the promoter. This increases the intensity of the extrinsic 
(repressor) noise in the low frequency regime. Moreover, 
the low-frequency fluctuations in the repressor binding do 
propagate through the different stages of gene expression. 
In particular, they lead to sharp bursts in the production 
of mRNA and protein. These bursts increase the noise 
intensity at the lower frequencies in the noise spectrum 
of mRNA and protein. And since the noise strength a^ is 
dominated by fluctuations in the low- frequency regime, 
spatial fluctuations ultimately strongly increase the noise 
in mRNA and protein concentration. 

Recently, experiments have been performed, in which 
the synthesis of individual mRNA transcripts |53| and 
individual protein molecules J53| could be detected. The 
systems in these studies were very similar to that stud- 
ied here: a gene under the control of a (Lac) repressor. 
These studies unambiguously demonstrated that tran- 
scription |53| and translation can occur in bursts [53J. 
Our simulation results show that spatial fluctuations of 
the repressor molecules might be responsible for this. In- 
deed, our results strongly suggest that spatial fluctua- 
tions are the dominant source of noise in gene expression 
in these systems. 

The spatial fluctuations due to diffusion of the repres- 
sor molecules could have significant implications for the 
functioning of gene regulatory networks. Under some 
conditions, it might be crucial that the protein number 
is not only low on average, but remains low at all times. 
For instance, if the protein itself functions as a transcrip- 
tion factor, it might by accident induce the expression of 
another gene, when, due to a fluctuation, its concentra- 
tion crosses a particular activation threshold. Thus, not 
all combinations of repressor copy number Nji and re- 
pressor backward rate fcbR that obey Eq. 1141 and thus 
have the same average repression strength, are neces- 
sarily equivalent in terms of function when diffusion is 
taken into account. If the fluctuations in the repressed 
state need to be small, then the cell could increase the 
number of repressors and decrease the binding affinity 
to the operator site, such that the repressor molecules 
stay bound to the DNA only briefly. Alternatively, the 
cell could minimize the effect of fluctuations by reducing 
the rate at which the open complex is formed by RNAP 
- our analysis shows that the process of open complex 
formation can act as a strong low-pass filter. 

The rapid rebindings observed in our simulations are a 
general phenomenon. We now address the question when 
the effect of spatial fluctuations due to diffusion can be 



13 



described by a well-stirred model in which the association 
and dissociation rates are renormalized. In the current 
problem, the rebinding time for a dissociated repressor 
is exceedingly short. As a consequence, the probability 
that a RNAP binds to the promoter during this time, is 
vanishingly small. This is precisely the reason that the 
effective dissociation rate is simply the bare dissociation 
rate divided by the number of rebindings (see Ea. ll8|l : the 
effective association rate is renormalized accordingly, be- 
cause the equilibrium constant should remain unchanged 
(see Ea. ll9f) . The success of the renormalized well-stirred 
model is thus a result of the strong separation of time 
scales - the time scale of repressor rebinding is well sep- 
arated from that of RNAP binding. 

The separation of time scales also makes it possible to 
account for the effect of spatial fluctuations by renormal- 
izing the association and dissociation rates in other cases. 
For instance, we have simulated a system in which repres- 
sion occurs in a cooperative manner. In this system, the 
repressor backward rate is smaller when two repressors 
are bound to the operator than when a single repressor 
is bound. However, when one of the two repressors disso- 
ciates, its rebinding time is so short that the probability 
for the other repressor to dissociate in the mean time, 
is negligible for reasonable values of cooperativity. As a 
result, the effect of spatial fluctuations can be described 
by a well-stirred model with properly renormalized reac- 
tion rates. We have also studied a system in which the 
expression of a gene is not under the control of a repres- 
sor, but rather under the control of an activator. Also 
in this system, diffusion of the transcription factors leads 
to an enhancement of noise in gene expression through a 
similar mechanism. 

Do these observations imply that the effect of spatial 
fluctuations can always be described by a well-stirred 
model? In the system studied here, the ligand (repressor) 
molecules bind to a single site. We expect that the ef- 
fect of spatial fluctuations becomes more intricate when 
the number of binding sites for a particular ligand in- 
creases - the binding of the ligand to the different sites 
will then exhibit correlations. This could be important 
when the ligand binds to receptors that occur in dense 
clusters, as in bacterial chemotaxis |5J, |53 and in the 
immune response |5q |. In gene regulatory networks this 
effect could also be significant. Recently, we have shown 
that in E. coli, pairs of co-regulated genes - genes that 
are controlled by a common transcription factor - tend 
to lie exceedingly close to each other on the genome [53 : 
their promoter regions are often separated by a distance 
smaller than a few hundred base pairs. It is conceivable 
that spatial fluctuations of the transcription factors in- 
troduce correlations between the noise in the expression 
of these pairs of co-regulated genes. This study also re- 
vealed that pairs of genes that regulate each other, often 
lie close together, again suggesting that the diffusive mo- 
tion of transcription factors could be important for the 
functioning of gene regulatory networks '53] ■ 

Even in the case of a single gene, the effect of spatial 



fluctuations is expected to be more subtle than that re- 
ported here. In this study, the operator is modeled as 
a spherical site. However, as mentioned in section III Al 
transcription factors are believed to find their operator 
site via a combination of free 3D diffusion and ID sliding 
along the DNA. While on length scales larger than the 
sliding distance this process is indeed essentially 3D dif- 
fusion, on length and time scales smaller than the sliding 
distance and sliding time, respectively, the dynamics is 
more complicated. We expect that sliding could have two 
important effects. First, it will increase the number of re- 
bindings - the probability that in ID a random walker 
returns to the origin is one, while in 3D there is a finite 
probability that it will escape and never return. Sec- 
ondly, sliding is expected to also increase the duration of 
the rebindings, especially when diffusion along the DNA 
is much slower than diffusion in the cytoplasm. It is thus 
conceivable that with sliding, the non-exponential relax- 
ation of the operator state, arising from the rebindings, 
shifts to lower frequencies (see Fig.QJ, where fluctuations 
in the operator state are not filtered out. In addition, it 
is possible that with sliding, RNAP and repressor com- 
pete for binding to the promoter on similar time scales, 
which would mean that the effective dissociation rate is 
no longer simply given by the bare dissociation rate di- 
vided by the number of rebindings. Indeed, under these 
conditions, the effect of spatial fluctuations is likely to 
become non-trivial, and describing it would probably re- 
quire a spatially resolved model. We leave this for future 
work. 



Finally, we address the question whether spatial fluctu- 
ations, and, more in particular, the rebindings, could be 
studied experimentally. Interestingly, recent biochemical 
data on the restriction enzyme EcoRV suggests that after 
an initial dissociation, 10-100 rebindings occur before the 
enzyme escapes into the bulk solution p35l l36l | , in good 
agreement with the average number of rebindings cal- 
culated in section IVTl However, in our gene expression 
model, the rebinding times are so short that it would 
seem difficult to probe the repressor rebindings directly 
in an experiment. In fact, reaction rates measured bio- 
chemically will probably already be corrected for accord- 
ing to Eos. 1161 and 1171 Sliding along the DNA, however, 
may extend the rebinding times to accessible experimen- 
tal time scales. Moreover, recent experiments suggest 
that the motion of proteins in the nucleoid might be sub- 
diffusive, which would increase the importance of the re- 
bindings 58J . Recently, magnetic tweezer experiments on 
a mechanically stretched, supercoiled, single DNA have 
made it possible to study the kinetics of the open complex 
formation and promoter clearance 26] . Performing these 
experiments on a promoter that is under the control of a 
repressor, seems a promising approach for studying the 
effect of spatial fluctuations due to the diffusive motion of 
transcription factors on the dynamics of gene expression. 
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The real and imaginary parts of the Fourier Transform 
are thus: 



^[X{lu)] = -ySkismcotk) (32) 

k 

^[X{lu)] = -Y^SkicosLotk), (33) 



The power spectrum of the time trace of the copy num- 
ber X{t) of a species X can be efficiently computed by 
exploiting the fact that in between the times tk the signal 
X{t) is constant. The Fourier Transform 5'x(w) oi X{t) 
is 

X{lu) = f X{t)e-"^'dt = J2 f^ Xke^^^'dt. (29) 



As X(t) is constant within every interval {tk-i,tk}, the 
integration can easily be performed: 






-iujtk-1' 



(30) 



Shifting up by one the index j in the second part of the 
sum, we obtain: 



X{u) = —y{Xk+i-Xk){t 

II,) ■^— ^ 



— iujtk 



(31) 



where we have defined 5k = ^fe+i — Xk- The Power 

~ 2 ~ 2 

spectrum Sy^{(jj) = 3ff [X(w)] + 3[X(w)] is thus given by 



S^{uj)={-y5kCos{ujtk)\ + -JI' 



I A; sm 



[ujtk) 



(34) 



The Fourier Transforms were computed at 10000 log- 
arithmically spaced angular frequencies starting from 
(^min = 10 • 27r/T, where T is the total length of the 
signal. Power spectra obtained according to Ea. 1341 were 
filtered with a box average over 20 neighboring points. 
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